import plotly.graph_objects as go
import plotly.figure_factory as ff
import plotly.express as px
import numpy as np
import pandas as pd

# plot force diagram
fig = go.Figure()
fig.add_shape(type='rect', x0=0, y0=1, x1=4, y1=3, fillcolor='rgba(133, 163, 201, 0.5)')
fig.add_scatter(x=np.linspace(-1,5, 20), y=np.ones(20)*3.5, line_color='black')
fig.add_scatter(x=np.linspace(-1,5, 20), y=np.ones(20)*0.5, line_color='black')
fig.add_annotation(x=-0.1, y=2, ax=-1, ay=2, xref='x', yref='y', axref='x', ayref='y',arrowhead=5, text=
                   r'$\require{color}{\color[rgb]{0.315209,0.728565,0.037706}p}$')
fig.add_annotation(x=4.1, y=2, ax=6, ay=2, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, text=
                   r'$\require{color} {\color[rgb]{0.315209,0.728565,0.037706}p} + \frac{\partial \
                   {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x} \delta x$')
fig.add_annotation(x=1, y=0.9, ax=3, ay=0.9, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, 
                   text=r'$\require{color}{\color[rgb]{0.990308,0.800015,0.121270}\tau}$')
fig.add_annotation(x=3.5, y=3.25, ax=0.5, ay=3.25, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, text=
                  r'$\require{color}{\color[rgb]{0.990308,0.800015,0.121270}\tau} + \frac{\partial \
                  {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y} \delta y$')
fig.add_scatter(x=[2], y=[3.65], mode='text', text=r'$\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}=0$')
fig.add_scatter(x=[2], y=[0.4], mode='text', text=r'$\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}=0$')
fig.add_annotation(x=-1.2, y=3.5, ax=-1.2, ay=0.4, xref='x', yref='y', axref='x', ayref='y', arrowhead=5)
fig.add_scatter(x=[-1.3], y=[2], mode='text', text=r'$H$')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

Because the velocity is constant in the x direction and depends only on the vertical position $F_{x}=0$ and because velocity is zero in the y direction $F_{y}=0$ \ The shear stress on the fluid is a function only of y and the pressure is a function only of x, therefore in the y-direction $\require{color}\frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial x} = \frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial y} = 0$ \ In order for the x components of shear and pressure to sum to zero they must both be constant along their respective axes \ By using the above force diagram we can get a relationship between $\require{color}\frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x}$ and $\require{color}\frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y}$: \ $$ \require{color}F_{x} = 0 = {\color[rgb]{0.315209,0.728565,0.037706}p}\delta y - \left({\color[rgb]{0.315209,0.728565,0.037706}p}+\frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x} \delta x\right)\delta y - {\color[rgb]{0.990308,0.800015,0.121270}\tau}\delta x + \left({\color[rgb]{0.990308,0.800015,0.121270}\tau} + \frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y} \delta y\right)\delta x \longrightarrow \frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x} = \frac{\partial {\color[rgb]{0.990308,0.800015,0.121270}\tau}}{\partial y} $$ Because pressure is only a function of x and shear stress is only a function of y we can say the partial derivatives from above can be converted to total derivatives: $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} = \frac{d{\color[rgb]{0.990308,0.800015,0.121270}\tau}}{dy}$ Furthermore we know the equation for shear stress: $\require{color}{\color[rgb]{0.990308,0.800015,0.121270}\tau} = {\color[rgb]{0.990448,0.502245,0.032881}\mu} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}}{dy}$ so we can get a relationship between velocity and the pressure gradient: $$ \require{color}\frac{d}{dy}\left({\color[rgb]{0.990448,0.502245,0.032881}\mu} \frac{d{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}}{dy}\right) = {\color[rgb]{0.990448,0.502245,0.032881}\mu} \frac{d^{2}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}}{dy^{2}} = \frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} $$ In Poiseuille flow we assume that the pressure gradient is constant but not zero, therefore we can solve the above equation for $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x}$ by assuming $\require{color}\frac{\partial {\color[rgb]{0.315209,0.728565,0.037706}p}}{\partial x}$ is some constant value thus the solution is in the form $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x} = \frac{1}{2{\color[rgb]{0.990448,0.502245,0.032881}\mu}}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}y^{2} + By + C$ where B and C are determined by the boundary conditions \ The boundary conditions are that at $\require{color}y=0\ {\color[rgb]{0.059472,0.501943,0.998465}v}=0$ and at $\require{color}y=H\ {\color[rgb]{0.059472,0.501943,0.998465}v}=0$ from these we get that $\require{color}{\color[rgb]{0.059472,0.501943,0.998465}v}_{x} = \frac{1}{2{\color[rgb]{0.990448,0.502245,0.032881}\mu}}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx}y(y - H)$

# define necessary variables 
H=10
dpdx = -1e-2
mu = 1e-3
x=np.zeros(11)
y=np.linspace(0, H, 11)
v_x= (1/(2*mu))*(dpdx)*y*(y-H)
v_y=np.zeros(11)

# create a vector plot showing the velocity vector of particles along y direction each with same x position 
fig = ff.create_quiver(x=x, y=y, u=v_x, v=v_y, arrow_scale=0.04)
fig.update_layout(xaxis=dict(range=[-1, 20]))
fig.add_scatter(x=np.linspace(-1,15, 20), y=np.ones(20)*10, line_color='black')
fig.add_scatter(x=np.linspace(-1,15, 20), y=np.zeros(20), line_color='black')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

# define how long it takes for the fastest particle to travel the length of the pipe
max_v = v_x[list(v_x).index(max(v_x))]
time = 20/max_v
distance = time*v_x

# using this time animate how the particles would idealy move in Poiseuille flow
fig = go.Figure()
frames = [go.Frame(data=[]) for k in range(50)]
for i, y_pos in enumerate(y[1:len(y)-1]):
    fig.add_scatter(x=np.linspace(0, 20, 20), y=np.ones(20)*y_pos, line_dash='dash', line_color='lightsteelblue')
    x_pos = np.linspace(0, distance[i+1], 50)
    for j,f in enumerate(frames):
        f.data += (go.Scatter(x=[x_pos[j]], y=[(np.ones(50)[j]*y_pos)], mode='markers', marker=dict(color='teal', size=10)),)

fig.update(frames=frames)

fig.update_layout(updatemenus = [dict(type='buttons', buttons=[dict(label='play', method='animate', 
                  args=[None, {'frame':{'duration':50}}])])], showlegend=False, xaxis=dict(range=[0,20]))

fig.add_scatter(x=np.linspace(0, 20, 20), y=np.ones(20)*H, line_color='black')
fig.add_scatter(x=np.linspace(0, 20, 20), y=np.zeros(20), line_color='black')
fig.update_layout(yaxis=dict(scaleanchor="x", scaleratio=1, showticklabels=False), xaxis=dict(showticklabels=False),
                  showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.show()

All of the above figures are based on the assumption that there is a negative pressure gradient, meaning the pressure gradient is favorable \ When the pressure gradient is greater than zero the pressure gradient is adverse If the flow were Couette+Poiseuille and the pressure gradient were adverse the flow would be pushed oppoiste to the direction of the moving plate near the bottom unmoving plate \ In just Poiseuille flow the the velocity profile reverses and the flow would move in the opposite direction to the initial velocity, which is in this case rightwards

# define an array of pressure gradients
dpdxs = np.linspace(-1e-2, 1e-2, 17)

# create a dictionary of all the y position and velocity data at each pressure gradient 
all_data = {}
for p in dpdxs:
    all_data[p] = {'y':[], 'v':[]}
    for y_val in y:
        v_xy = ((1/(2*mu))*(p)*y_val*(y_val-H))
        all_data[p]['y'].append(y_val)
        all_data[p]['v'].append(v_xy)

# create a plot showing how velocity profile changes as the pressure gradient changes 
fig = go.Figure()
fig.add_scatter(x=np.linspace(-13, 13, 20), y=np.ones(20)*H, line_color='black')
fig.add_scatter(x=np.linspace(-13, 13, 20), y=np.zeros(20), line_color='black')

frames = []
for k in range(0,len(dpdxs)):
    figaux = ff.create_quiver(x=x, y=all_data[dpdxs[k]]['y'], u=all_data[dpdxs[k]]['v'], v = v_y, arrow_scale=0.04, 
                              line_color='steelblue')
    frames.append(go.Frame(data=[figaux.data[0]]))

fig.update_layout(updatemenus=[dict(type='buttons', buttons=[dict(label='Play', method='animate',
                  args=[None, dict(frame=dict(duration=500, redraw=True))])])], xaxis=dict(range=[-13,13], showticklabels=
                  False), yaxis=dict(showticklabels=False), showlegend=False, plot_bgcolor='rgba(0, 0, 0, 0)')
fig.update(frames=frames)

for k in range(len(fig.frames)):
    fig.frames[k]['layout'].update(title_text=f'Pressure Gradient {round(dpdxs[k],3)}')

fig.add_scatter(x=np.linspace(-13, 13, 20), y=np.ones(20)*H, line_color='black')
fig.add_annotation(x=5, y=-1, ax=-5, ay=-1, xref='x', yref='y', axref='x', ayref='y', arrowhead=5, text = r'$+x$')

fig.show()

In the above animation you can see that when $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} > 0$ the flow is moving to the left opposing the initial rightwards velocity of the flow \ When the flow instead has $\require{color}\frac{d{\color[rgb]{0.315209,0.728565,0.037706}p}}{dx} < 0$ the flow moves rightwards in the same direction as the initial velocity

# It would be good to add a horizontal axis showing what direction positive x is.